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We determine the limiting dynamics of a fermionic condensate following a sudden perturbation 
for various initial conditions. We demonstrate that possible initial states of the condensate fall into 
two classes. In the first case, the order parameter asymptotes to a constant value. The approach to 
a constant is oscillatory with an inverse square root decay. This happens, e.g., when the strength 
of pairing is abruptly changed while the system is in the paired ground state and more generally 
for any noneqmUbrium state that is in the same class as the ground state. In the second case, 
the order parameter exhibits persistent oscillations with several frequencies. This is realized for 
noncquilibrium states that belong to the same class as excited stationary states. Our classification 
of initial states extends the concept of excitation spectrum to nonequilibrium regime and allows one 
to predict the evolution without solving equations of motion. 



The response of a fermionic condensate to fast exter- 
nal perturbations presents a long-standing problenijli- 
H^ . The main difficulty is to describe the time evolu- 
tion in the non-adiabatic regime when a nonequilibrium 
state of the condensate is created on a time scale shorter 
than the Cooper instability time ta = I/Aq, where Ag 
is the equilibrium BCS gap. In this case the evolution of 
the system cannot be described in terms of a quasiparti- 
cle spectrum or a single time-dependent order parameter 
A{t) 13]. One has to account for the dynamics of in- 
dividual Cooper pairs, making it a complex many-body 
problem. 

The non-adiabatic regime can be accessed experimen- 
tally in ultra-cold Fermi gases, where the strength of pair- 
ing between fermions can be rapidly changed 14). Non- 
adiabatic measurements can be also performed in quan- 
tum circuits utilizing superconducting qubits (nanoscale 
superconductors). Here nonequilibrium conditions can 
be generated by fast voltage pulses on a time scale com- 
parable to TA |l5|. 

Here we consider a BCS condensate that is out of equi- 
librium at t = and study its time evolution for i > 0. 
Given the state of the system at < = 0, we predict the dy- 
namics with no need for actually solving equations of mo- 
tion. We show that physically meaningful initial states 
fall into two broad categories. In the first scenario, |A(t)| 
asymptotes to a constant value A^o < Aq. The approach 
to Aqo is oscillatory with a decay, 



1 



cos(2Aooi 



(1) 



where the constants a and (j) depend on details of the ini- 
tial state. This scenario is realized, for example, when the 
pairing strength is abruptly changed, while the system is 
in the paired ground state. In the second scenario, |A(t)| 



oscillates persistently with several incommensurate fre- 
quencies. The number of frequencies as well as the limit- 
ing dynamics of individual pairs can be inferred directly 
from the initial state. 

We propose a topological classification of initial states, 
which extends the concept of excitation spectrum to the 
nonequilibrium regime. If a state is in the same class as 
the paired ground state, Eq. applies. Other states 
are topologically distinct, in which case persistent oscil- 
lations occur. 

Our approach explains differences between previous 
studies of condensate dynamics. A linear analysis around 
the BCS ground state yieldsJ3| damped oscillations with 
a frequency 2Ao. Eq. |^ generalizes this result to 
the nonlinear case and a wide range of initial condi- 
tions. An oscillatory decay following a chan ge i n the 
coupling strength was observed numerically [slTlH . We 
will see that this is due to the fact that initial states of 



Refs. are in the same class as the BCS ground 

state. Undamped periodic oscillations of I A(t)| have been 
found in Refs. ^IMIHl' They were also seen numerically 
for initial states close to a normal state0- In contrast, 
Ref. |T3 also starts from the normal state , bu t obtains 
a saturation to A^o = Ao/2. It turns out[l6j that this 
occurs if the initial state is a paired state with a small 
seed gap A^ <^ An- Quasiperiodic oscillations of the 
order parameter [9l [lOj can also be realized (see below). 

In the non-dissipative regime, one can use the BCS 
model to describe the dynamics of the condensate. Here 
we are interested in the thermodynamic limit. Then, the 
BCS mean-field is valid as long as the order parameter 
is nonzero. There are several equivalent ways to derive 
mean-field evolution equations. Using Anderson's pseu- 
dospin representationjlj, one can describe the mean-field 
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evolution by a classical spin HamiltonianQ 
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where ej are the single-particle energies and Sj—Sji: 
is^. Dynamical variables Sj are vectors of fixed length, 
\sj\ ^ 1/2. The BCS order parameter is A{t) =^ ~ 
iAy — g'^jSj. Mean-field equations of motion sj — 
{Hbcs, Sj} follow from Hamiltonian (O and the angular 
momentum Poisson brackets, {sj, s^} = —SjkSj etc., 



Sj = hj X Sj hj = {—2Ax, —'2Ay, 2ej) 



(3) 



The general non-stationary BCS problem can now be for- 
mulated as follows: given an initial distribution of spins, 
determine its evolution. 

Solution of Eq. ^ for each of the dynamical variables 
Sj{t) and A(<) is a superposition of harmonic oscillations 
with certain frequencies. In the thermodynamic limit 
the frequency spectrum in general consists of continuum 
and discrete parts. Consider, e.g., the absolute value of 
the order parameter |A(i)|. Suppose the discrete part 
contains k incommensurate frequencies. Let their contri- 
bution to |A(t)| be Fkit), 
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where D is the ultraviolet cutoff (Debye frequency in 
conventional BCS). The integral in Eq. |0J vanishes as 
t — > oo for any well behaved A{uj). Hence, | A(i)| Fk{t) 
for t oo, i.e. |A(t)| oscillates with k frequencies for 
fc > 1 and \A{t)\ ^ Aoo if fc = 0. 

Given the initial state, the structure of the frequency 
spectrum can be determined using integrability of BCS 
dynamics 0, 01 . Frequencies of an integrable system de- 
pend only on its integrals of motion|l7|. which can be 
evaluated at t — 0. Frequencies can thus be determined 
without solving equations of motion. 

Here we use the method of Ref. . It is convenient 
to introduce the following vector function (Lax vector) 
of an auxiliary (spectral) parameter u: 
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where z is a unit vector along the z-axis. The square 
of the Lax vector is conserved by the evolution. The 
frequency spectrum is related to branch cuts of w{u) = 
yJlJ^u). Note that the numerator of (u) is a polyno- 
mial of degree 2n, where n is the total number of lev- 
els tj. Since L^(m) > 0, all 2n roots come in complex 
conjugate pairs. For finite n, all roots are typically dis- 
tinct leading to n isolated cuts connecting pairs of con- 
jugate roots. This situation is described by the general 
solution^ for finite n. One can also have a situation 



when 2(n — m) roots are real and therefore double de- 
generate. This leaves m branch cuts corresponding to 
TO pairs of complex conjugate roots. The dynamics can 
now be described in terms of only m < n effective spins 
governed by Hamiltonian (j^J with to spins and m new 
effective energy levels. These m-spin solutions contain 
only m incommensurate frequencies. One frequency cor- 
responds to a uniform rotation of all spins around the 
z-axis. Thus, |A(t)| contains m — 1 frequencies |l(lj |. 

In the thermodynamic limit some roots of (u) merge 
into continuous lines and give rise to the continuum part 
of the spectrum, while isolated pairs of roots correspond 
to the discrete part. Thus, the number of discrete fre- 
quencies in \ A{t)\ is the number of isolated pairs of roots 
less one, fc = m — 1. At large times |A(i)| exhibits persis- 
tent oscillations with fc frequencies |l9| and is described 
by an TO-spin solution J^. The number fc is a topological 
property of the initial state. It is the number of handles 
on the Riemann surface of the function y'L^u). 

Discrete part of the frequency spectrum turns out to be 
related to discontinuities of the spin distribution s(e) as a 
function of e. To see this, consider first stationary states. 
There are two types of such states. The BCS ground state 
and excited states with a constant A 7^ are obtained by 
aligning each spin Sj self-consistently along its effective 
magnetic field hj . These states can be termed anomalous 
stationary states. Choosing the x-axis so that A is real, 
we obtain 
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where ej = 1 if the spin is antiparallel to the field and 
Cj = —1 otherwise. The self-consistency condition A — 
9 ^j yields the BCS gap equation. The state with all 
Cj = 1 is the BCS ground state. The state ek = —1 and 
ej^k = 1 is a state with a single excited pair(l| of energy 
2\/£fe + ^0- Using Eqs. (j3El) and the gap equation, we 
derive 



Ls(7i) — —(Ax + Uz)is(w), 
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We see that L^(u) = {u'^ + A'^ )Ll{u) has a pair of isolated 
roots at M = ±zA. All other roots are determined by the 
equation Ls{u) — and are double degenerate. 

First, consider the ground state, Cj = 1. Since Ls(m) — > 
±00 as w — > EjiO for each j, all roots of (u) are real and 
located between consecutive ej. In the thermodynamic 
limit they merge into a line from ~D to D as shown in 
Fig. Note that the existence of a double real root 
between ej and ej+i relies on ej ~ ej+i. Further, Eq. © 



implies that when e,- 



the components of spins 



s(e) are continuous at e = ej. Now let Cj 



Bj+l. 



In 
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this case the real root between cj and Cj+i can disap- 
pear. Thus, discontinuities (jumps) in the spin distribu- 
tion generate isolated complex roots (see Fig. ^ for an 
example). Since spins far from the Fermi level are not 
flipped, the total number of jumps is even. In general, 
for 2p jumps Ls{u) can have up to p pairs of isolated 
roots. 

One can study equations of motion ||2J) linearized 
around anomalous stationary states. Setting 5sj = 
Aj-e*'^*'*, we solve for the normal modes. The eigenvalues 
turn out to be LOj = 2a/?x| + A^, where Uj are the roots 
of L^(u). For the ground state Xj — ej up to finite size 
corrections (cf. Ref. 0). 

Next, consider few examples of initial states far from 
equilibrium. Let the system be in an anomalous station- 
ary state for t <0. Suppose att — the coupling changes 
abruptly from g' to g. Using Eq. (jSJ, one can show that 
the change in g results in a smooth deformation of the 
root distribution. Lines of roots deform into lines. On 
the other hand, doubly degenerate roots become non- 
degenerate. A state that had p pairs of degenerate iso- 
lated roots in addition to a pair of roots ±i A now has 
m = 2p+l pairs of non-degenerate roots, i.e. 2p+l cuts 
of yX^(u). As shown above, in this case |A(t)| exhibits 
persistent oscillations with k ~ m ^ I ~ 2p frequencies. 
This behavior is illustrated in Fig. 

Let the initial state be the ground state with coupling 
g' . Then, the line of double real roots splits into two 
complex conjugate lines (Fig. ^p). There is only one 
pair of isolated roots as in the ground state. Therefore, 
k = and |A(i)| asymptotes to a constant Aao at i S> ta, 
as illustrated in Fig. According to Eq. (0), at large 
times spin Sj rotates in a constant magnetic field hj — 
(2Aoo, 0, 2ej) with a frequency uj{ej) — 2\/e| -I- A^. Us- 
ing this, one derives Eq. (^). The 1/^/t decay law is set 
by the square root singularity in the spectral densitvflij. 
Note that, even though |A(i)| asymptotes to a constant, 
the final state of the system is non-stationary. In the final 
state each spin precesses with its own frequency [l^. 

There is another type of stationary states - normal 
states. In these states each spin is aligned along the z- 
axis, Sj = Zj/2 = ±1/2. The Fermi state is Zj = — sgne^ 
(levels below the Fermi energy are occupied, above - 
empty). States with other Zj correspond to particle-hole 
excitations of the Fermi gas. For example, a state Zfe — 
sgn Cfe < has a pair of fermions removed from the level 
Efc. The Lax vector for normal states is L„(w) — L„{u)z, 
where 

L„iu) = -^-+J2-^ (9) 

All roots of L^(w) — L\{u) are thus doubly degenerate. 
Note the absence of a branch cut of ^/\]^u) connecting 
the points u = ±iA. Further analysis is similar to that 
for anomalous states. The Fermi state has a single jump 



in the Sz{e) at the Fermi level. This results in a pair of 
complex conjugate isolated roots (Fig.l^t-), which can be 
determined from the equation i„(u) = 0. In the thermo- 
dynamic limit, we obtain u = ±iAo/2. The rest of the 
roots are real and form a line from —D to D. For a gen- 
eral normal state with 2p+l iumDs[2o| in Sz(e), one can 
have up to p -|- 1 complex conjugate pairs of roots. Each 
root is doubly degenerate. Linearizing Eqs. © around 
a normal state, we obtain normal frequencies loj = 2uj, 
where Ln{uj) = 0. In particular, the Fermi state has a 
single unstable mode that corresponds to uj = ±iAo/2. 
This mode grows as e^"* indicating the pairing instability 
of the Fermi state 0,|23|. Remaining frequencies are real 
and correspond to the precession of spins at their natural 
frequencies, tj^ — 2efc up to finite size corrections. 

Let the system be in or close to a normal state at 
t — 0. First, we consider the Fermi state. Since, within 
mean-field, normal states are unstable equilibria, a small 
perturbation is needed to start off the dynamics. One can 
start e.g. from a non-stationary spin distribution close to 
the Fermi statejji. A typical deviation splits all double 
degenerate roots as illustrated in Fig.|2|3. Real roots split 
into two complex conjugate lines close to the real axis. 
Their contribution to | A(t)| is small in the deviation from 
the Fermi state and decaying. Degenerate complex roots 
at u = ±iAo/2 split into m = 2 isolated cuts close to each 
other. Since fc = m— l = lin Eq. Q, |A(t)| will exhibit 
undamped periodic oscillations, as shown in Fig|3;. Its 
functional form is described by a 2-spin solution ,18.j. The 
roots close to the real axis indicate that the initial spin 
distribution does not match the 2-spin solution exactly 
(see the discussion of m-spin solutions above). 

The dynamics in normal states can also be triggered 
by quantum fluctuations. In Ref. this is modelled by 
adding a small external fleld — 2Aqfx to Eq. ||2J). This 
violates conservation of L^(u). We have d'L'^{u)/dt = 
2/S.QFLy/ g. Being applied for a short time t* , the exter- 
nal fleld drives the system out of the normal state. Treat- 
ing the evolution of L^(m) perturbatively, we find that the 
new positions of the roots are determined by the equa- 
tion Ln{u) — ±2i/S.QFt* / g. Degenerate complex roots 
split along the real axis into two cuts by 2AQpt*d/g, 
where d = {cj+i — ej) is the level spacing. The resulting 
root configuration shown in Fig. [Sja contains two cuts. 
Thus, k — 1 and |A(<;)| exhibits periodic oscillations, see 
Fig.Et. 

We see that for initial conditions close to the Fermi 
state a periodic solution is "dynamically selected" ■ 
Other initial conditions "select" damped oscillations or 
multi- frequency undamped oscillations. In more conven- 
tional terms, this corresponds to a basic fact that the 
evolution of an integrable system depends on initial con- 
ditions. All diferent behaviors are captured by the gen- 
eral solution^]. Here we systematically classified possible 
initial states and specialized the general solution to each 
type of initial conditions, i.e. we developed "selection 
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rules" for the BCS dynamics. Which behavior is reahzed 
in a particular experimental setup depends on the initial 
state of the condensate. In this respect, the periodic so- 
lution is somewhat special; if we start from the ground 
state with a small Ai„ ^ Aq, the order parameter |A(t)| 
asymptotes to a constant value Aqq, see Eq. (Q. In the 
thermo dyna mic limit, the damping disappears only when 
A„, = 00. 

Excited normal states have several jumps in the spin 
distribution and can therefore display oscillations with 
more than one frequency. Consider, e.g., a state where 
pairs in energy interval from — Eq to —tb have been re- 
moved. The initial spin distribution is Sz(e) = — sgn[e(e-|- 
£«)(£ + £&)]• The 2p+l = 3 jumpsHOl result mp+l = 2 
pairs of isolated double degenerate complex roots. As a 
small perturbation splits these roots into m = 2p -I- 2 = 4 
cuts, |A(i)| oscillates with k = m— l = 2p-\-l — i 
incommensurate frequencies, see Fig. 0] 

In conclusion, we have shown how the non-stationary 
BCS dynamics can be predicted from the initial state of 
the condensate. We classified initial states by their in- 
tegrals of motion - roots of L^(w). For states with a 
root diagram as in the paired ground state, the order 
parameter |A(i)| displays damped oscillations described 
by Eq. Other states are of the same type as ex- 

cited stationary states of the BCS Hamiltonian. In these 
cases |A(t)| oscillates persistently with few incommensu- 
rate frequencies. The number of frequencies is related 
to the number of jumps in the pseudospin distribution 
of the corresponding stationary state. For the situation 
most relevant to the experiments on cold fermions - an 
abrupt change of the coupling - we predict damped os- 
cillations with \/\/t decay. 
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FIG. 1: a) Roots Uj of L'^(m) for the BCS ground state of n = 100 spins with equilibrium gap Aq. Axes: Im(itj)/Ag and 
Re(«j)/d, where d is the level spacing. Insets: spin distribution s^(e) and s^(e) (s''(e) = 0). The absence of discontinuities 
in the spin distribution (2p = 0) implies a single pair of isolated roots at iiAg (red circles connected by a dashed line). 
Remaining 2n — 2 = 198 roots (solid green circles) are real, doubly degenerate, and located between consecutive energy levels 
ej (ticks on the real axis). 

b) At t = the coupling constant is abruptly increased so that the corresponding ground state gap is Aq = 2.4Aq, 
while the spin configuration remains the same as in part a). The line of double real roots deforms into two complex 
conjugate lines. The Riemann surface of •y^L^(it) has m = 2p + 1 = 1 isolated cuts (red circles connected by a dashed line). 
This nonequilibrium state has the same number of isolated roots as the BCS ground state and, therefore, is in the same 
topological class. The frequency spectrum of |A(t)| is thus continuous: k — m— 1 = 2p = 0, i.e. no discrete frequencies in Eq. 



c) Time evolution of |A(t)| after the change Aq — > Aq. In the absence of discrete frequencies, the order parameter 
|A(t)| asymptotes to a constant value Aoo (see Eq. (Q). 
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FIG. 2: a) Roots u-, of L^(ii) for the Fermi state of n = 100 spins with ground state gap Aq. Axes: Im(itj)/Ao and Re(iij)/d, 
where d is the level spacing. Insets: spin distribution s^(e) — — sgne/2 (s'^(e) — s^(e) = 0). There is 2p + 1 = 1 discontinuity 
(a jump in s^(e) at the Fermi level) in the spin distribution resulting in p + 1 = 1 pair of imaginary doubly degenerate roots 
Uj = ±iAo/2 (blue crosses). Remaining 2n ~ 4 — 196 roots (solid green circles) are real, doubly degenerate and located 
between consecutive energy levels ej (ticks on the real axis). Since all spins are along the z-axis, this state is stationary (see 
Eq. The normal frequencies are Ljj — 2uj. The mode corresponding to Uj — ±iAo/2 is unstable, it grows as e^"'. 

b) Initially the system is in the Fermi state. A small "external magnetic field" — 2Aqfx = —5.4 x 10~^Aox is added 
to Eq. Q for a time t* = I/Aq. The p + 1 = 1 pair of double roots Uj = ±jAo/2 splits into m = 2p + 2 = 2 pairs of isolated 
roots (red circles connected by dashed lines); the line of real roots splits into two complex conjugate lines. This nonequilibrium 
state has the same number of isolated roots as the Fermi state and, therefore, belongs to the same topological class. The 
Riemann surface of ^yT^Ju) has m = 2p + 2 = 2 isolated cuts. The frequency spectrum of |A(t)| has k = m~l — 2p + l — 1 
discrete frequency in Eq. I^J- 

c) Time evolution of |A(t)| for the initial state described in part b). Since there is one discrete frequency, |A(f)| 
exhibits periodic oscillations. Its functional form is described by an (m = 2)-spin solution (see the text below Eq. 
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FIG. 3: a) Roots Uj of L^(u) of n = 100 spins in an excited stationary state. This anomalous state (see the text) has been 
obtained from the ground state with the gap Ag by flipping spins in energy interval (— 0.37Ao,0). It has a gap A = 0.4Aq. 
Axes: Im(uj)/Ao and Re(«j)/d, where d is the level spacing. Insets: spin distribution s^(e) and s^(e) (s^(e) = 0). Spin flips 
result in 2p = 2 discontinuities in the spin distribution. Therefore, there is a p = 1 pair of isolated doubly degenerate roots 
(blue crosses) in addition to a pair of roots Uj = ±iA (red circles connected by a dashed line). The remaining 2n — 6 = 194 
roots (solid green circles) are real, doubly degenerate, and located between consecutive energy levels ej (ticks on the real axis). 

b) At t = the coupling constant is abruptly increased so that the corresponding ground state gap is Ao = 1.55Aq, 
while the spin conflguration remains the same as in part a). The p = 1 pair of isolated double roots splits into 2p = 2 pairs 
of isolated roots. Together with a pair of roots coming from Uj = ±iA, there are m = 2p + 1 = 3 pairs of isolated roots (red 
circles connected by dashed lines). The line of double real roots deforms into two complex conjugate lines. The Riemann 
surface of yJ'L'^{u) has m = 2p + 1 = 3 isolated cuts. This nonequilibrium state has the same number of isolated roots as the 
excited state in part a) and therefore is in the same topological class. The frequency spectrum of |A(f)| has fc = m — l = 2p = 2 
discrete frequencies. 



c) Time evolution of |A(t)| after the change Aq Aq. Since there are k = m 
Eq. 01 the order parameter |A(t)| displays oscillations with two basic frequencies, 
(m = 3)-spin solution (see the text below Eq. 0). 
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FIG. 4: a) Roots itj of L^(ii) for an excited stationary state of n = 100 spins with ground state gap Aq. Axes: Im(tij)/Ao 
and Re(uj)/d, where d is the level spacing. This normal state (see the text) has been obtained from the Fermi state by flipping 
spins in energy interval (— 0.5Ao, — Aq). Insets: spin distribution s^(e) (s^(e) = s''(e) = 0). There are 2p+ 1 = 3 discontinuities 
(three jumps in s^(e)) in the spin distribution resulting in p + 1 = 2 pairs of isolated double roots Uj — fj,\,2 ± 271,2 (blue 
crosses). Remaining 2n — 8 = 192 roots (solid green circles) are real, doubly degenerate and located between consecutive 
energy levels e-, (ticks on the real axis). Since all spins are along the z-axis, this state is stationary (see Eq. The 
normal frequencies are ujj = 2uj. There are two unstable modes corresponding to Uj = /ii,2±i7i,2, which grow as e^''^* and e^''^*. 

b) Initially the system is in the stationary state described in part a). A small "external magnetic field" 
— 2Aqfx — —5.4 X 10~^Aox is added to Eq. ^ for a time t* — I/Aq. The p + 1 = 2 pairs of degenerate roots 
Uj — /ii,2 ± 171,2 split into m = 2p + 2 = 4 pairs of isolated roots (red circles connected by a dashed line), while the line 
of real roots splits into two complex conjugate lines. This nonequilibrium state has the same number of isolated roots as 
the stationary state in part a) and, therefore, belongs to the same topological class. The Riemann surface of y'L^u) has 
m = 2p + 2 = 4 isolated cuts. Insets: Fourier spectra of ReA(t) and ImA(t) display m = 2p + 2 = 4 basic frequencies. One 
frequency (vo) corresponds to a uniform rotation of all spins around the z-axis and cancels out in |A(t)|. Thus, the frequency 
spectrum of | A(t)| has fc = m — l = 2p+l = 3 discrete frequencies. 

c) Time evolution of |A(t)| for the initial state described in part b). Since A: = m — l = 2p-+-l = 3in Eq. 
|A(t)| exhibits quasiperiodic oscillations with three basic frequencies (see insets in Fig. ^t)- Us functional form is described 
by an (m — 4)-spin solution (see the text below Eq. 0). 



